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ABSTRACT 

Flow vorticity is a fundamental property of turbulent convection in rotating systems. Solar supergranules exhibit a preferred sense 
of rotation, which depends on the hemisphere. This is due to the Coriolis force acting on the diverging horizontal flows. We aim to 
spatially resolve the vertical flow vorticity of the average supergranule at different latitudes, both for outflow and inflow regions. To 
measure the vertical vorticity, we use two independent techniques: time-distance helioseismology (TD) and local correlation tracking 
of granules in intensity images (LCT) using data from the Helioseismic and Magnetic Imager (HMI) onboard the Solar Dynamics 
Observatory (SDO). Both maps are corrected for center-to-limb systematic errors. We find that 8-h TD and LCT maps of vertical 
vorticity are highly correlated at large spatial scales. Associated with the average supergranule outflow, we find tangential (vortical) 
flows that reach about 10ms _1 in the clockwise direction at 40° latitude. In average inflow regions, the tangential flow reaches the 
same magnitude, but in the anti-clockwise direction. These tangential velocities are much smaller than the radial (diverging) flow 
component (300 m s -1 for the average outflow and 200 m s -1 for the average inflow). The results for TD and LCT as measured from 
HMI are in excellent agreement for latitudes between -60° and 60°. From HMI LCT, we measure the vorticity peak of the average 
supergranule to have a full width at half maximum of about 13 Mm for outflows and 8 Mm for inflows. This is larger than the spatial 
resolution of the LCT measurements (about 3 Mm). On the other hand, the vorticity peak in outflows is about half the value measured 
at inflows (e.g. 4 x 10 -6 s -1 clockwise compared to 8 x 10 -6 s -1 anti-clockwise at 40° latititude). Results from the Michelson Doppler 
Imager (MDI) onboard the Solar and Heliospheric Observatory (SOHO) obtained in 2010 are biased compared to the HMI/SDO 
results for the same period. 
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1. Introduction 


Duvall & Gizon|p000| and|Gizon & Duval~I| ( 


supergranules (see Rieutord & Rincon[2010 


2003] ) revealed that 
for a review) pos¬ 


sess a statistically preferred sense of rotation that depends on so¬ 
lar latitude. In the northern hemisphere, supergranules tend to ro¬ 
tate clockwise, in the southern hemisphere anti-clockwise. This 
is due to the Coriolis force acting on the divergent horizontal 
flows of supergranules. For supergranulation (lifetime > 1 day), 
the Coriolis number is close to unity (see |Gizon et al.|20l0| ». As 
a consequence, the vorticity induced by the Coriolis force should 
be measurable by averaging the vorticity of many realizations of 
supergranules at a particular latitude. 

For single realizations, |Attie et aL| ( |20Q9| ) detected strong 
vortices associated with supergranular inflow regions by apply¬ 
ing a technique called balltracking. Komm et al. |2007| pre¬ 
sented maps of vortical flows in quiet Sun convection using 
helioseismic ring-diagram analysis. With the same technique, 
Hindman et al. ( |2009| ) resolved the circular flow component as¬ 


sociated with inflows into active regions. The spatial structure of 
such vortical flows has not yet been studied though for many re¬ 
alizations. Knowledge of the flow structure of the average super¬ 
granule will help constrain models and simulations of turbulent 
convection that take into account rotation. 

Here, we aim to spatially resolve the vertical component of 
flow vorticity associated with the average supergranule. We in¬ 
vestigate both outflows from supergranule centers and inflows 


into the supergranular network. To measure the flow divergence 
and vorticity, we use two independent techniques: time-distance 
helioseismology (TD) and local correlation tracking (L CT) o f 
granules. We use the TD method from Langfellne r et al.| ( [2014| ), 
where a measurement geometry that is particularly sensitive to 
the vertical component of flow vorticity was defined. 


1.1. Time-distance helioseismology 

Time-distance helioseismology makes use of waves travelling 
through the Sun ( [Duvall et al.||[1993| ). A wave travelling from 
the surface point iq through the solar interior to another surface 
point r 2 is sensitive to local physical conditions (e. g. the wave 
speed or density). A flow in the direction r 2 - iq will increase 
the wave speed, thus reducing the travel time r + from iq to r 2 . A 
flow in the opposite direction will result in a longer travel time. 
The travel time is measured from the temporal cross-covariance, 
labeled C, of the observable 0 obtained at the points iq and r 2 : 

h N 

C{ ri,r 2 , t) = J F 0( r i> ti)(p(r 2 , t, + t), (1) 

i=-N 

where h t is the temporal cadence, T = 2(N + 1 )h t is the observa¬ 
tion time, and /)• = (- N, -N + 1,..., N)h t are the times when 0 is 
sampled. Typically, the observable 0 is the Doppler line-of-sight 
velocity component. 
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Fig. 1 . Measurement geometries for divergence- and vorticity-sensitive 
travel times and LCT velocities, a) The divergence-sensitive travel time 
r 01 is obtained by measuring the travel-time difference between a central 
point r and a surrounding annulus of radius A. b) The vorticity-sensitive 
travel time r ac is obtained by measuring the travel-time differences r dlff 
between adjacent points along a regular polygon surrounding r. The 
polygon consists of n points and has edges of length A. The points are 
situated on a circle of radius R = A/[2 sin(^/w)T The travel time r ac is 
the average over the components r dlff (see Eq. We obtain circula¬ 
tion velocities v ac from multiplying r ac by a calibration factor, c) From 
LCT, we obtain the horizontal velocity components v* and v y . Given a 
reference point r, the LCT velocities can be expressed in 2D polar coor¬ 
dinates (r, 6) by the outward pointing radial velocity component v r and 
the anti-clockwise pointing tangential velocity component v t . d) Using 
LCT velocities, we approximate v ac by averaging the tangential velocity 
component v t over the annulus shaded in green. The annulus is defined 
by its radius R and half-width s. We choose R = A for n - 6 and 
s = 2 Mm. 


The travel time can be obtained from C by fitting a wavelet 
( Duvall et al.||1997| ) or by comparison with a reference cross¬ 
covariance and application of an appropriate weight function 
( Gizon & Birch|2004| . To distinguish the flow signal in the travel 
time from other perturbations (e. g. local sound speed changes), 
we use the travel-time difference 


T xWt (r 1 ,r 2 ) = T + (r 1 ,r 2 )-r + (r 2 ,r 1 ). 


( 2 ) 


Travel times that are especially sensitive to the horizontal 
flow divergence can be obtained by replacing r 2 with an annu¬ 
lus around ri (see Fig.[lJ). Averaging (/> over the annulus yields 
the “outward-inward” travel time r 01 (Duvall et al. 1996). To ob¬ 


tain travel times that measure the vertical component of the flow 
vorticity, we average r dlff components along a closed contour in 
anti-clockwise direction ( [Langfellner et al.|2014| ). We choose the 
contour to be a regular polygon with n points and edge length A 
in order to approximate an annulus (see Fig. The mean over 
the r dlff components gives the vorticity-sensitive r ac travel time: 


r ac (r, A, n) := - V T dl V„r, + i), 
n 

i=i 

where we use the notation r n+ i = iq. 


( 3 ) 


1.2. Local correlation tracking 

LCT measures how structures in solar images are advected due 
to background flows. For the tracer, it is common to use solar 
granulation observed in photospheric intensity images ( |Novem-| 
|ber & Simon|1988| ). The general procedure is as follows. Pairs of 
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images are selected that observe the same granules but are sepa¬ 
rated by a time At. This time separation must be small compared 
to the lifetime of granules, i.e. At <$c 10 min. To obtain spa¬ 
tially resolved velocity maps, an output map grid is defined. For 
each grid point, subsets of the intensity images that are centered 
around the grid point are selected by applying a spatial window 
and multiplied by a Gaussian with a full width at half maximum 
(FWHM) of typically a few megameters. The subsets are then 
cross-correlated in the two spatial image dimensions x and y. The 
peak position (Ax, Ay) of the cross-correlation yields the spatial 
shift. Since the measured shift is usually only a small fraction of 
a pixel, it must be obtained using an appropriate fitting proce¬ 
dure. Finally, the velocity components in x and y directions are 
given by v* = Ax/At and v-y = Ay/At. 

The LCT method has prov en valuable to measure flo w pat- 
terns in the Sun. For instance, Bran dt et~aL| ( |1988| ) and |Simon| 
et al.|(|1989|) ob served single vortex flows at granulation scale. 


Hathaway et al.| ( [20T3| ) detected giant convection cells with LCT 


of supergranules in Doppler velocity images. For a comparison 
of different LCT techniques, see [Welsch et al.| ( [2007] ). 


2. Observations and data processing 

The basis for our measurements of wave travel times and flow 
velocities from local correlation tracking are two independent 
observables. We use Doppler velocity images for the TD and 
intensity images for the LCT. Both observables are measured 
for the full solar disk by the Helioseismic and Magnetic Imager 
(HMI) onboard the Solar Dynamics Observatory (SDO) ( |Schou| 
et al.|2012| ) and are available for the same periods of time. This 
allows a direct comparison of the two methods for looking at 
“the same Sun” but utilizing independent data. 

We used 112 days of both SDO/HMI Dopplergrams and in¬ 
tensity images in the period from 1 May through 28 August 
2010. Patches of approximate size 180x 180 Mm 2 were selected 
that are centered at solar latitudes from -60° to 60° in steps of 
20°. They were tracked for 24 h each at a rate consistent with 
the solar rotation rate from [Snodgrass| ( ^984| at the center of the 
patch. The data cubes cross the central meridian approximately 
at half the tracking time. They were remapped using Postel’s 
projection with a spatial sampling of 0.5 arcsec px -1 (0.348 Mm 
px -1 ). The temporal cadence is 45 s. We divided each data cube 
into three 8 h datasets. The x direction of the remapped images 
points to the west, the y direction points to the north. 

For further comparison, we also used Dopplergrams from the 
Michelson Doppler Imager (MDI) onboard the Solar and Helio¬ 
spheric Observatory (SOHO) spacecraft ( |Scherrer et al.||1995| ). 
We chose 59 days of images taken in the MDI full-disk mode 
that overlap in time with the HMI data (8 May through 11 July 
2010). We tracked and remapped the MDI Dopplergrams in the 
same manner as for HMI, albeit with a coarser spatial sampling 
of 2.0 arcsec px -1 (1.4 Mm px -1 ). 

2.1. Flow velocity maps from local correlation tracking 

For the LCT, we used our own code with the HMI photospheric 
intensity images as input. Our code is similar to the FLCT code 
by |Fisher & Welsch| ( [2008[ ), but uses another procedure to mea¬ 
sure the peak positions of the cross-correlation (described later 
in this section). 

We removed the temporal mean image for every dataset 
and chose an output grid with a sampling of 2.5 arcsec px -1 
(1.7 Mm px -1 ), thus five times coarser than the input images. 
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Fig. 2. Comparison of line-of-sight velocity from two different data products at 40° solar latitude on 6 June 2010. a) HMI Dopplergram averaged 
over 8 h. The map has been convolved with a Gaussian of cr/ V2 ~ 1.4 Mm and subsampled to match the coarser LCT resolution. The mean 
over the map and a linear function in the x direction (parameters determined by a least-squares fit) have been subtracted, b) LCT map from HMI 
intensity images, averaged over 8 h. The line-of-sight velocity component was computed from the v* and components. For v y , the mean over the 
map and a linear function in y direction have been subtracted, c) Scatter plot of the two maps. The Pearson correlation coefficient is 0.94. The red 
line shows the direction of largest scatter and crosses the origin. It is a best-fit line in the sense that it minimizes the sum of squared distances of 
the points perpendicular to the line (cf. |Pearson|1901| . This is different from linear regression, where no error in the x coordinate is assumed and 
only the sum of squared distances in the y coordinate is minimized. The slope of the red line is 1.08; the error in the direction of lowest scatter is 
33.9 m s _1 . 


The size of the image subsets used for the cross-correlation is 
adapted to the width of the Gaussian the subsets are multiplied 
by. We chose cr = 2 Mm for the Gaussian and a diameter of 
4<x = 8 Mm for the subsets both in x and y directions. The sub¬ 
sets are separated in time by At = 45 s (the cadence), which is 
sufficiently small compared to the granules’ evolution timescale. 
We averaged the cross-correlations over the whole 8 h dataset. 

To measure the peak position (Ax, Ay) of the cross¬ 
correlation, we calculated (separately for x and y directions) the 
parameters of a parabola matching the cross-correlation at the 
maximum and the adjacent pixels. To improve the estimate of the 
peak position, we translated the cross-correlation by (-Ax, -Ay) 
using Fourier interpolation and iterated the parabolic fit. We re¬ 
peated this procedure four times in total. The measured shifts 
converge quickly, the maximum additional shift in a fifth iter¬ 
ation is of the order 10 -5 px at 60° latitude (corresponding to 
0.2 ms" 1 or less), the root mean square of the additional veloc¬ 
ity shift is less than 0.02 m s _1 . The measured peak position is 
the sum of the shifts measured in each step. 

In Fig. |2j we compare the line-of-sight component vlos of 
the LCT velocity with the velocity of an average Dopplergram, 
obtained by averaging Dopplergrams over the same time-period 
as the LCT maps (8 hours). Both images are for the same re¬ 
gion at 40° latitude around the central meridian. At this latitude, 
the average Dopplergrams are dominated by the horizontal flows 
that can be measured with LCT, but systematic effects like fore¬ 
shortening are weak (see Appendix |C.1| ). We convolved the av¬ 
erage Dopplergram with a Gaussian of width cr/ V2 « 1.4 Mm 
(FWHM roughly 3 Mm). This resembles the convolution of the 
intensity maps prior to computing the correlation of image sub¬ 
sets in the LCT. The chosen width maximizes the correlation co¬ 
efficient between the average Dopplergram and the LCT image. 
In addition, we interpolated the average Dopplergram onto the 
coarser LCT grid and at each pixel subtracted the mean velocity 
over the map. To remove the residual rotation signal, we further 
subtracted a linear gradient in the x direction that we obtained 
from a least-squares fit of v x averaged over y. The LCT line- 
of-sight velocity component has been computed from v* and v y . 
The Vy map showed a linear gradient in the y direction leading 


to an average velocity difference of about -200 m s -1 between 
the bottom and the top of the map. This gradient is presumably 
due to the “shrinking Sun” effect, which has been discussed in 
Lisle & Toomre ( 2004), albeit for LCT of Dopplergrams (a short 
description is also given in Appendix |C.1| ). The gradient (and the 
mean over the Vy map) was removed before computing vlos • 
The processed vlos maps from direct Doppler data and LCT 
agree well (correlation coefficient 0.94). The scatter plot shows 
that on average the velocity values from LCT are slightly larger 
(by a factor 1.08) than from the Dopplergrams. This is differ¬ 
ent from what other authors have reported. De Rosa & Toomre 


(1998) measured a slope of 0.89 and later (De Rosa & Toomre 


|2004| ) 0.69 using their LCT co de and SOHO/MD I Doppler 
grams. [Rieutord et aL| ( |2001| ) and |Verma et al.| ( [2013| ) found that 
LCT underestimates the real velocities in convection simulations 
(however at smaller spatial scales than we are studying). 


2.2. Travel-time maps for horizontal divergence and vertical 
vorticity 

As input for the travel-time measurements, we used the HMI 
Dopplergrams. In Fourier space, we filtered the 8 h datasets to 
select either the f-mode or pi-mode ridge. The filters consist 
of a raised cosine function with a plateau region in frequency 
around the ridge maximum for every wavenumber k. Addition¬ 
ally, power for kR Q < 300 and kR Q > 2 600 (f modes) and for 
kR Q < 180 and kR Q > 1800 (pi modes) respectively is dis¬ 
carded. The symbol R e denotes the solar radius. The filter details 
are given in Appendix [A| 

We computed the cross-correlation C from each filtered 
Doppler dataset in temporal Fourier space using 


C(r u r 2 ,oj) = hu<f>* T (r u <o)(f> T (r 2 ,a)), (4) 

where oj denotes the angular frequency, h w = In/T (with T = 
8 h) is the frequency resolution, and (j) T is the temporal Fourier 
transform of the filtered dataset (multiplied by the temporal win¬ 
dow function). 

For each dataset, we measured travel times r 01 with an an¬ 
nulus radius of 10 Mm and T ac with the parameters A = 10 Mm 
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and n = 6 (regular hexagon). We rotated the hexagon structure 
successively three times by an angle of 15° to obtain four r ac 
measurements for the same dataset that are only weakly corre¬ 
lated (see |Langfellner ^t al.||2014| for details). Averaging over 
these measurements yields a higher signal-to-noise ratio (S/N). 
We used the linearized travel-time definition by Gizon & Birch 
( |2004| ) and a sliding reference cross-covariance that we obtained 
by averaging C over the entire map. 


3. Comparison of horizontal divergence and 
vertical vorticity from TD and LCT 


We now want to compare the measurements of horizontal diver¬ 
gence and vertical vorticity from TD and LCT for HMI. For TD, 
we use r 01 and r ac as divergence- and vorticity-sensitive quan¬ 
tities. For LCT, we can directly compute div/* = d x v x + d y v y 
and oj z = d x v y - d y v x from the v* and v } , maps. To compute the 
derivatives of the LCT velocities, we apply Savitzky-Golay fil¬ 
ters ( jSavitzky & Golay 1964| for a polynomial of degree three 
and a window length of 15 pixels (about 5 Mm, with a FWHM 
of about 3 Mm of the smoothing kernel). The Savitzky-Golay 
filters smooth out variations in the derivatives on spatial scales 
below the LCT resolution. 


In the case of vorticity, we can also attempt a more direct 
comparison of TD and LCT. Consider the horizontal velocity 
field in 2D polar coordinates around r (see Fig. [1]:). Instead of 
using v x and v y for LCT, we can study the velocity component in 
the radial (divergent) direction, v r , and the component in the tan¬ 
gential (anti-clockwise) direction, v t . The travel time r ac essen¬ 
tially measures v t averaged over the closed contour. The travel¬ 
time r ac is built up of point-to-point components r dlff that capture 
the flow component that is parallel to the line connecting the two 
measurement points. The velocity magnitude that corresponds 
to the travel time r dlff can roughly be estimated by calibration 
measurements using a uniform flow (Appendix |B|. We use this 
calibration to convert r ac travel times into flow velocities and 
call the result v ac . Since convective flows are highly turbulent, a 
conversion factor obtained from uniform flows has to be treated 
with caution though. Additionally, the conversion factor is sen¬ 
sitive to the details of the ridge filter (Appendix |C. 3 1 ). Also note 
that because no inversion is applied, the velocities v ac represent 
an average over a depth range given by travel-time sensitivity 
kernels. For f modes, the range is from the surface to a depth 
of about 2 Mm, with a maximum of sensitivity near the surface, 
and for pi modes from the surface to roughly 3 Mm, with one 
maximum near the surface and another one at a depth of about 
2 Mm (see, e.g., |Birch & Gizon|2007| . With LCT, we approxi¬ 
mate v ac by averaging v t over an hard-edge annulus with radius 
R = 10 Mm and half-width s = 2 Mm (see Fig.[l]l). The annulus 
width roughly corresponds to the width of travel-time sensitivity 
kernels (see, e.g.,[Jackiewicz et al. 2007] ). 


For the divergence-sensitive measurements, such a compari¬ 
son is not possible without an inversion of the r 01 maps. There¬ 
fore, we stick to comparing TD r 01 and LCT div/* in the follow¬ 
ing. 


3.1. Spatial power spectra of horizontal divergence and 
vertical vorticity 

From the TD r 01 and r ac maps as well as the LCT di.v/* and co z 
maps, we calculated the spatial power spectra and averaged them 
over azimuth. The result for HMI is shown in Fig. [3] Note that 
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we rescaled the amplitude of the LCT power in order to show it 
together with the travel-time power. 

For the divergence, the TD and LCT powers show a simi¬ 
lar behavior at larger scales (except for kR Q = 25, which corre¬ 
sponds to the map size). However, all three curves peak at differ¬ 
ent scales - f modes at kR Q = 120, pi modes at kR Q = 100 and 
LCT at kR Q = 150. The comparison with the curves for the TD 
noise model ( [Gizon & Birch|2004| shows that the highest S/N 
for the TD r 01 occurs at supergranulation scale, with pi modes 
probing slightly larger scales than f modes. For LCT, no noise 
model is available that we know of. Thus it remains unclear if the 
peak of the power coincides with the peak of the S/N. For small 
scales ( JcRq larger than 300) the LCT power vanishes quickly, 
whereas the TD power reaches a noise plateau (f at kR e = 500, 
pi at kR Q = 300). 

In the case of vorticity, the curves for TD and LCT look 
similar at large scales, albeit the power for LCT oj z drops more 
quickly toward larger scales than for TD r ac . Compared to the 
divergence case, the peak positions are slightly shifted toward 
larger scales. However, the comparison with the TD noise model 
reveals that the S/N does not have a peak at supergranulation 
scale but continues to increase toward larger scales (cf. Langfell- 
|ner et al.||2014| ). At mid scales, the LCT power drops off only 
slowly, whereas the TD power quickly reaches the noise level (f 
at kR Q = 250, pi at kR Q = 200). It is not clear if the considerably 
larger power of LCT co z at mid scales (150 < kR Q < 500) is due 
to real flows or noise. At smaller scales, both TD curves behave 
more erratically. This happens, however, in a regime of almost 
pure noise. LCT power drops off quickly beyond kR Q = 400. 


3.2. Maps of horizontal divergence and vertical vorticity 

For comparing maps of horizontal divergence and vertical vor¬ 
ticity, one point to consider is the different spatial sampling for 
TD and LCT maps. To correct for this, we interpolate the ve¬ 
locity maps derived from LCT onto the (finer) travel-time grid. 
In order to compare the maps on different spatial scales, we ap¬ 
ply different band-pass filters to the individual maps in Fourier 
space. The individual filters are centered around kR Q values of 
50 through 400 in steps of 50. Each filter is one in a plateau re¬ 
gion of width 50, centered around these values. Adjacent to both 
sides of the plateau are raised cosine flanks that make the fil¬ 
ter smoothly reach zero within a kR Q range of 50. Additionally, 
we employ a high-pass filter for kR Q > 400. From all maps, we 
subtract the respective mean map over all 336 datasets prior to 
filtering. 

Example 8 h maps for r 01 and v ac from f-mode travel times 
as well as di.v/* and co z from LCT are depicted in Fig. [4] The 
maps are filtered around kR Q = 100. Note that for the sake of 
an easier comparison, we plotted -r 01 rather than r 01 . For the 
flow divergence, all three maps are highly correlated. The aver¬ 
age correlation coefficients over all 336 maps are 0.96 between 
LCT div/* and -r 01 for f modes and 0.92 between LCT div/* and 
-r 01 for pi modes. 

In the case of flow vorticity, the agreement of the LCT and 
TD maps is weaker than for the divergence. The average cor¬ 
relation coefficient over all 336 maps is 0.68 between LCT oj z 
and f-mode v ac and 0.51 between LCT oj z and pi-mode v ac (not 
shown). For comparing LCT v ac instead of co z with TD v ac , the 
correlation coefficients are noticably higher (0.75 for f modes 
and 0.57 for pi modes). The flow magnitudes are roughly com¬ 
parable. 

Table Q] shows the correlation coefficients between LCT and 
TD averaged over all datasets for all filters and including pi 
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Fig. 3. Power spectra (averaged over az¬ 
imuth of wavevector k and 336 datasets) of 
TD travel-time maps and LCT d\v h and oj z 
maps computed from HMI Dopplergrams 
and intensity images at the solar equator 
near disk center (distributed symmetrically 
between 7° east and west of the central 
meridian), a) Divergence-sensitive travel 
times r 01 for f modes and pi modes as well 
as LCT div/j. b) Vorticity-sensitive travel 
times r ac for f modes and pi modes as well 
as LCT (x) z . The amplitudes of LCT di\ h 
and (x) z power were rescaled to match the 
range of the travel-time power. The dashed 
lines represent noise models for the f and 
Pi modes based on |Gizon & Birch| ( |2004| ). 
The thickness of the lines denotes the lcr 
error. 


Fig. 4. Comparison of TD and LCT maps 
at the equator. All maps are based on 8 h 
of HMI images (intensity and Doppler ve¬ 
locity) taken on 1 May 2010 and have been 
band-pass-filtered around kR Q = 100. The 
colorbar limits are set to the maximum ab¬ 
solute value of the corresponding map and 
symmetrized around zero. Left column: 
Divergence-sensitive travel times r 01 for f 
modes as well as LCT horizontal diver¬ 
gence div/j. Right column: Circulation ve¬ 
locities v ac for TD as well as LCT ver¬ 
tical vorticity oj z . The TD r ac maps have 
been converted into velocity maps by point- 
wise multiplication with a constant factor 
-5.62 m s -2 (f modes) and -11.1 ms -2 (pi 
modes), see Appendix |B| for details. 


modes. The error in the correlation coefficients is less than 0.01. 
Note that the edges (12 Mm) were removed from the maps be¬ 
fore the correlation coefficients were computed. For the flow di¬ 
vergence, the correlation coefficients are almost constantly high 
for smaller kR Q values. In the range kR Q = 300 - 400, the 
correlation coefficient between LCT div/* and -r 01 for f modes 
rapidly decreases from 0.78 to 0.31. For LCT and pi modes, 
the correlation coefficient decreases from 0.83 to 0.15 from 
kR Q = 200 - 350. For the high-pass filters, the LCT and TD 
maps are completely uncorrelated. 

In the case of vorticity, the correlation decreases rapidly for 
both f and pi modes at kR Q = 200. Again, the LCT and TD 
maps are uncorrelated for large kR Q . The correlation coefficients 
for LCT v ac are significantly higher than for LCT co z . 

The dependence of the correlation coefficients on spatial 
scale conceptually agrees well with the power spectra in Fig. [3] 


There is a high correlation on large scales where the observed 
TD travel-time power clearly exceeds the power of the TD noise 
model. On the other hand, the very low correlation on smaller 
scales reflects that the power of TD observations and noise 
model are almost equal. 

Qualitatively, the correlation coefficients are comparable 
with the value 0.89 from [De Rosa et aL] ( )2000| who obtained 
travel-time and LCT velocity maps from SOHO/MDI Doppler¬ 
grams and smoothed the divergence maps by convolving with a 
Gaussian with FWHM 6.2 Mm. 


4. Net vortical flows in the average supergranule 

The major goal of this paper is to spatially resolve the vorticity 
of the average supergranule at different solar latitudes. In the 
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Fig. 5. Correction for geometrical center-to-limb systematics in maps of circulation velocity v ac for LCT and TD using HMI observations. This 
example shows the correction of v ac maps at 60° latitude, a) Patch 1 (raw data) is located at 60° latitude at the central meridian. Patch 2 (systematics 
data) is situated at the equator, 60° west of the central meridian. We correct the raw data by subtracting patch 2 from patch 1 (after rotating patch 
2 by 90° anti-clockwise). Note that we also use measurements of the systematics from the east for the correction (rotated 90° clockwise). This is 
omitted in the sketch for clarity, b) Maps of v ac for LCT and f-mode TD. The three columns of panels show (from left to right) the measured raw 
data at 60° latitude, the measured systematics from 60° west and east (rotated and averaged), and the corrected data. 


Table 1. Correlation between LCT maps and TD travel-time maps de¬ 
rived from HMI intensity images and Dopplergrams. 




Correlation coeff. between LCT and TD 

Modes 

kRo 

LCT di \h 

LCT 0J Z 

LCT v ac 

(TD) 


TD -r oi 

TD -r ac 

TD —r ac 

f 

50 

0.93 

0.70 

0.77 


100 

0.96 

0.68 

0.75 


150 

0.96 

0.63 

0.68 


200 

0.94 

0.53 

0.57 


250 

0.89 

0.31 

0.30 


300 

0.78 

-0.01 

0.14 


350 

0.58 

-0.09 

0.23 


400 

0.31 

-0.06 

0.23 


>400 

0.02 

0.00 

0.03 

Pi 

50 

0.90 

0.53 

0.59 


100 

0.92 

0.51 

0.57 


150 

0.89 

0.44 

0.50 


200 

0.83 

0.33 

0.38 


250 

0.66 

0.13 

0.19 


300 

0.36 

-0.06 

0.11 


350 

0.15 

-0.05 

0.11 


400 

0.04 

-0.02 

0.08 


>400 

0.00 

0.00 

0.01 


following, we describe the averaging process and show average 
divergence and vorticity maps. 

4.1. Obtaining maps of the average supergranule 

To construct the average supergranule, we started by identifying 
the location of supergranule outflows and inflows in f-mode t 01 
maps from HMI and MDI. We smoothed the maps by remov¬ 
ing power for kR Q > 300 and applied an image segmentation 
algorithm (cf. [Hirzberger et al.|2008| ). The coordinates of the in¬ 


dividual supergranules were used to align maps of various data 
products. For each identified position, we translated a copy of 
the map to move the corresponding supergranule to the map cen¬ 
ter. These translated maps were then averaged. At each latitude, 
we averaged over roughly 3 000 supergranules in total for HMI 
(1500 supergranules for MDI). Supergranules closer than 8 Mm 
to the map edges were discarded. 

We produced maps for the average supergranule outflow and 
inflow from r 01 and r ac travel-time maps as well as LCT v x and 
v_y maps. Prior to the averaging process, the LCT maps were spa¬ 
tially interpolated onto the (finer) travel-time grid. For all maps, 
we subtracted the mean map over all 336 HMI datasets (177 
datasets in the case of MDI). This removes signal that does not 
(or only slowly) change with time, including differential rotation. 
Additionally, we removed power for kR e > 300 by applying a 
low-pass filter in Fourier space. 

The resulting average r ac maps were converted into v ac maps. 
From the LCT v* and v y maps for the average supergranule, 
we computed div/*, oj z , and v ac . The Savitzky-Golay filters that 
we employed to compute the spatial derivatives smooth out step 
artefacts from the image alignment process, yet preserve the sig¬ 
nal down to the resolution limit of the LCT. 

We corrected the v ac and co z maps for geometrical center-to- 
limb systematics (unless stated otherwise). We measured these 
effects using HMI and MDI observations west and east of the 
disk center, at relative longitudes corresponding to the latitudes 
of the regular observations. The idea is that any difference (be¬ 
yond the noise background) between maps at disk center and 
a location west or east from disk center is due to geometrical 
center-to-limb systematics. Such systematics only depend on the 
distance to the disk center. Therefore, our raw measurements of 
v ac and (x) z that we obtained north and south of the equator should 
be affected by the systematics in the same way as measurements 
west and east of the disk center. We corrected the raw data by 
subtracting the v ac and oj z maps west and east of the disk center. 
This approach is analogous to |Zhao et al.| ( |2013| who used the 
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method to correct measurements of the meridional circulation. 
Figure [5]illustrates the correction process for v ac maps at 60° lat¬ 
itude. The correction is particularly important for LCT at high 
latitudes. We note that the measured center-to-limb systematics 
at lower latitudes (up to 40° north and south) are much weaker 
and only lead to a mild correction of the v ac and oj z maps. A fur¬ 
ther discussion of the center-to-limb systematics can be found in 
Appendix [CJ] 


4.2. Latitudinal dependence of the vertical vorticity in outflow 
regions 

Figure [6] shows the circulation velocity v ac in the average su¬ 
pergranule outflow region for LCT and f-mode TD for latitudes 
from -60° to 60°, in steps of 20°. For comparison, the left col¬ 
umn shows the horizontal divergence div/* from LCT. At all lati¬ 
tudes, there is a peak of positive divergence at the origin. All di¬ 
vergence peaks are surrounded by rings of negative divergence. 
This reflects that on average every supergranule outflow region 
is isotropically surrounded by inflow regions. The strength of the 
divergence peak slightly decreases toward higher latitudes. Fur¬ 
thermore, the divergence peaks are slightly shifted equatorward 
at high latitudes (by about 0.7 Mm at ±60°). These effects are 
presumably due to center-to-limb systematics. 

The v ac maps (center and right columns) show negative peaks 
(clockwise motion) in the northern hemisphere and positive 
peaks (anti-clockwise motion) in the southern hemisphere. The 
peaks are surrounded by rings of opposite sign, as for the di¬ 
vergence maps. There is a remarkable agreement between LCT 
and TD in both shape and strength of the peak structures. At 
the solar equator, no peak and ring structures are visible. Note 
that the LCT and TD v ac maps at the equator are still correlated 
though. This shows that the “noise” background is due to real 
flows rather than measurement noise that is dependent on the 
technique. 

To study the latitudinal dependence of the observed and cor¬ 
rected signal in more detail, we plot in Fig. [7^ the peak velocity 
v ac from Fig.[6j including pi-mode TD, as a function of solar lat¬ 
itude (lines). The peak velocity shows an overall decrease from 
south to north, with a zero-crossing at the equator. The curves 
are antisymmetric with respect to the origin. The peak velocities 
have similar values at a given latitude, with f-mode velocities 
appearing slightly stronger than LCT and pi-mode velocities (in 
this order). The highest velocities are slightly above 10 m s -1 . 
Figure [7]} shows the peak magnitude in maps of the vertical vor¬ 
ticity oj z , as measured from LCT. The overall appearance is simi¬ 
lar to the circulation velocities v ac . The highest absolute vorticity 
value is about 5 x 10 -6 s -1 . 

Figure [8^ shows cuts through y = 0 for the maps of LCT and 
TD v ac (including pi modes) at 40° latitude. We choose this lat¬ 
itude because the S/N in the v ac and co z peaks is high compared 
to other latitudes, whereas the measurements are only mildly af¬ 
fected by center-to-limb systematics. The velocity magnitudes 
and shapes of the curves are comparable for the three cases. For 
the LCT and f-mode curves, an asymmetry in the west-east di¬ 
rection is visible. This means that the ring structures surround¬ 
ing the peaks in the v ac maps are stronger in the west than in the 
east. The FWHM is about 13 Mm in all cases. The peaks are 
very slightly shifted eastwards. However, this east shift does not 
appear to be a general feature at all latitudes. Mostly, the shifts 
are consistent with random fluctuations. Partly, the shifts might 
also be due to other effects, for instance an incomplete removal 
of center-to-limb systematics. 


LCTdiv h A LCT v ac TD v ac (f modes) 


60° 


40° 


20 ° 


0 ° 


- 20 ° 


-40° 


-60° 


-40-20 0 20 40 -40-20 0 20 40 -40-20 0 20 40 

x [Mm] 


-200 0 200 -15 -10 -5 0 5 10 15 

div h [icrv 1 ] v ac [ms -1 ] 

Fig. 6. Maps of horizontal divergence div /2 and circulation velocity v ac 
for the average supergranule outflow regions at various solar latitudes 
(after correction for geometrical center-to-limb systematics). The maps 
have been derived from HMI intensity images and Dopplergrams. The 
horizontal divergence div /2 has been computed from LCT v x and hor¬ 
izontal velocity maps. The LCT v ac maps were obtained by averag¬ 
ing the tangential velocity component v t over an annulus with radius 
10 Mm and half-width 2 Mm to resemble the r ac measurement geom- 
etry (Fig. [ljl). For TD, the v ac maps are based on r ac travel-time maps 
(Fig.[]J>) that have been computed from Dopplergrams. 


40 

20 

0 

-20 

-40 






For comparison, the FWHM of the t 01 peak structure is about 
13 Mm for pi modes, compared to about 11 Mm for f-mode r 01 . 
The horizontal divergence div/ 2 from LCT at 40° latitude peaks 
at about 170 x 10 -6 s -1 with a FWHM of about 10 Mm. 
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Fig. 7. Peak v ac and oj z values for the 
average supergranule at different solar 
latitudes, a) Circulation velocity v ac for 
LCT, f modes and pi modes, b) Verti¬ 
cal component of flow vorticity co z ob¬ 
tained from LCT. Solid lines are for 
the average supergranule outflow re¬ 
gion, dashed lines for the average super¬ 
granule inflow region. At 0° latitude, the 
values at the map center are shown in¬ 
stead of the peak values. The errorbars 
have been computed from dividing the 
336 datasets into eight parts and mea¬ 
suring the variance of v ac and oj z at the 
peak positions over the eight parts. 



Fig. 8. a) Cuts through the maps of 
the circulation velocity v ac for the aver¬ 
age supergranule outflow region at 40° 
latitude (shown in Fig. [6j, at y = 0. 
The TD and LCT maps were derived 
from HMI Dopplergrams and intensity 
images. The thin lines denote estimates 
of the variability of the data as obtained 
from dividing the 336 datasets into eight 
parts. The lcr level is shown, b) As a), 
but for the average supergranule inflow 
region at 40° latitude. 


From the v ac peak velocities, we can estimate the average 
vorticity (oj z )a over the circular area A of radius R = 10 Mm that 
is enclosed by the r ac measurement contour (see Fig. [TJ. The 
average vorticity is given by 


(^z)a 


r 

A 


2v ac 

~~R~ 


( 5 ) 


where T is the flow circulation along the r ac measurement con¬ 
tour that we approximated with T ~ 2nRv ac . By taking the 
v ac peak values, we obtain (qj z )a ~ -2.4 x 10 -6 s -1 for the f 
modes, (oj z )a ~ -1.6 x 10 -6 s -1 for the pi modes, and (oj z )a ~ 
-2.Ox 10 -6 s -1 for LCT. Thus the average vorticity in the circular 
region is roughly half the peak vorticity at 40° latitude. 


4.3. Inflow regions 

So far, we have discussed vortical flows around supergranule 
outflow centers. It is interesting though to compare the magni¬ 
tude and profile of these flows with the average inflow regions, 
which have a different geometrical structure (connected network 
instead of isolated cells). Analogue to Fig. [6] for the outflows, 
Fig. [9] shows maps of di\ h and v ac around the average supergran¬ 
ule inflow center. As for the outflows, the v ac maps from TD 
and LCT agree very well at all analyzed latitudes. The peaks 
in the v ac maps have the opposite sign compared to the outflows. 
This indicates that flows are preferentially in the clockwise (anti¬ 
clockwise) direction in the average supergranular outflow region 
and anti-clockwise (clockwise) in the average inflow region in 
the northern (southern) hemisphere. Cuts through y = 0 of the 
v ac maps at 40° latitude are shown in Fig. J8J). The v ac curves 
have the same shape as the corresponding curves for the aver¬ 
age outflow center (with a FWHM of 14 to 16 Mm) but the peak 
flow magnitude is reduced and the sign is switched. As for the 
outflows, the ring structures are stronger on the west side than 
on the east side. 
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The horizontal flow divergence div/* in the average inflow is 
similar to the average outflow (about the same FWHM) but with 
reversed signs and reduced magnitude. The peak divergence is 
about -120 x 10 -6 s _1 at 40° latitude with a FWHM of about 
10 Mm. As for the outflows, there is a systematic decrease in 
peak magnitude and a slight equatorward shift of the diV/* peak 
at high latitudes. 

The latitude dependence of the v ac peak values for the aver¬ 
age supergranule inflow region (dashed lines in Fig. [7^) is almost 
mirror-symmetric to the outflow regions. The values are slightly 
smaller though compared to the average outflow, with a ratio in¬ 
flow/outflow of -0.87 ± 0.03 for the f modes, -0.85 ± 0.06 for 
the pi modes, and -0.72 ± 0.05 for the LCT v ac . In the case of oj z 
(Fig. [ 7 ])), on the other hand, the ratio between the average inflow 
and outflow center is -1.8 ± 0.2. 

From the peak values of v ac , we can estimate the average 
vorticity (co z )a over the circular area A of radius R = 10 Mm 
in the same way as for the outflow regions. We obtain (cj z )a ~ 
1.8 x 10 -6 s -1 for the f modes, (oj z )a ~ 1.2 x 10 -6 s -1 for the pi 
modes, and (o) z )a ~ 1.4 x 10 -6 s -1 for LCT. The peak vorticity at 
40° latitude is therefore larger by a factor of about five compared 
to the vorticity averaged over the circular area. 


4.4. Dependence of the vertical vorticity on horizontal 
divergence 

The detection of net tangential flows in the average supergranule 
raises the question of how much the magnitudes of v ac and oj z 
depend on the selection of supergranules. As a test, we sort the 
identified supergranules at 40° latitude from HMI with respect 
to their divergence strength, as measured by the peak f-mode r 01 
of each supergranule. The sorted supergranules are assigned to 
four bins, which each contain roughly the same number of su¬ 
pergranules. The boundaries of the bins for f-mode r 01 are about 
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Fig. 9. Same as Fig. [6] but for the average supergranule inflow regions. 


-96.9, -53.8, -42.1, -31.5, and -16.0 s for the outflows and 
67.5, 38.3, 32.1, 26.0, and 11.1 s for the inflows. Note that a 
simple scatter plot would be very noisy because the v ac and cj z 
maps are dominated by turbulence. 

For each bin, we computed the peak TD r 01 and r ac as well as 
LCT di \h and co z in the same way as for all identified supergran¬ 
ules that we discussed in the previous sections, but without the 
correction for center-to-limb systematics. In Fig.|T0[i, we plot the 
peak r ac as a function of the peak r 01 from f modes and pi modes 
both for outflows and inflows. The magnitude of r ac clearly in¬ 
creases with r 01 . The ratio r ac /r 01 is roughly constant. Only the 
f-mode bin for the weakest inflows deviates substantially from 
this behaviour. Figure \TUjp shows the peak co z versus the peak 


di \h from LCT. In this case, the relationship is less clear, con¬ 
sidering the large vertical errorbars. A constant ratio oj z /divi t is, 
at least, by eye consistent with the measurements. However, for 
outflows oj z might also be constant. Note that the fit lines for the 
travel times in Fig.fTOfc have almost the same slopes for outflows 
and inflows, whereas in the case of LCT div/* and co z the slope for 
the inflows is much steeper than for the outflows. This is consis¬ 
tent with Fig. [ 7 ] where cj z was shown to be twice as strong in the 
inflows as in the outflows, whereas the velocities v ac are of sim¬ 
ilar mag nitud e (n ot ju st for TD, but also for LCT). As discussed 
in Sects. [42| and |4.3[ the velocities v ac do not directly measure 
the vorticity at a given position but rather a spatial average. 

In general, we can conclude that a selection bias in favour 
of stronger or weaker supergranules probably does not affect the 
measured ratio of vertical vorticity to horizontal divergence. 


4.5. Comparison of SDO/HMI and SOHO/MDI 

While the results for the average supergranule have been ob¬ 
tained using different methods (TD and LCT) and image types 
(Dopplergrams and intensity images), they are all based on the 
same instrument - HMI. It is thus useful to compare the HMI re¬ 
sults to v ac maps that have been measured from independent MDI 
data. Since MDI cannot sufficiently resolve granules at higher 
latitudes to successfully perform LCT, however, we only discuss 
TD. 

In contrast to HMI, the correction for geometric center-to- 
limb systematics is not sufficient for MDI. For example, for f- 
mode TD at 40° latitude, the central peak structure appears elon¬ 
gated (see Fig. |C.l| in the appendix). Nevertheless, the v ac values 
at the origin are remarkably similar for HMI and MDI. At the 
average outflow, we measure (-11.1 ± 0.4) m s -1 (HMI) versus 
(-10.1 ± 0.8) m s _1 (MDI) for f modes and (-7.7 ± 0.3) m s _1 
compared to (-6.5 ± 0.8) m s -1 for pi modes. 

For inflows, the MDI v ac maps compare to HMI in the same 
manner, with MDI being slightly weaker than HMI. The flow 
magnitudes for HMI and MDI at the origin after correction are 
(8.9 ± 0.4) m s _1 versus (7.2 ± 1.0) m s _1 for f modes and (6.1 ± 
0.3) m s -1 compared to (4.2 ± 0.5) m s -1 for pi modes. Note that 
the noise background is stronger for MDI. This is, however, not 
surprising, since only about half the number of Dopplergrams 
(compared to HMI) have been used to produce these maps. 

The latitude dependence of v ac at the origin for MDI is quali¬ 
tatively comparable with HMI (see Fig. [IT]). We measure a zero¬ 
crossing and sign change of v ac at the equator, both for the av¬ 
erage supergranule outflow and inflow regions. However, the v ac 
magnitudes are systematically smaller for MDI. This difference 
increases farther away from the equator. It is especially dramatic 
for f modes at ±60° latitude. Whereas v ac reaches values between 
10 and 12 ms' 1 at these latitudes in HMI, for MDI the velocity 
magnitudes are below 5 ms" 1 . This is probably connected to the 
lower spatial resolution of MDI, which results in a larger impact 
of geometrical foreshortening effects at high latitudes compared 
to HMI. 

While MDI clearly performs much worse than HMI, the 
agreement with HMI at the origin gives reason to believe that 
MDI v ac measurements are still useful. This would be especially 
interesting for long-term studies of the solar cycle dependence 
since continuous data reaching back to 1996 could be used. 


5. Differences between outflow and inflow regions 

The differences between the average supergranule outflow and 
inflow regions as measured from LCT in HMI data are summa- 
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Fig. 10. Vorticity as a function of divergence 
for the average supergranule at 40° latitude 
from HMI data, a) Vorticity travel times r ac 
versus divergence travel times r 01 for TD f 
modes and pi modes. The upper-left quadrant 
shows the values for the average outflows, the 
lower-right quadrant for the average inflows. 
The solid and the dashed lines show least- 
squares fits of a linear function through the 
origin for outflows and inflows, respectively. 
The errorbars were obtained from dividing the 
336 datasets into eight parts, b) As a) but the 
peak oj z versus the peak 6\\ h from LCT is 
shown. Note that the quadrants depicting out¬ 
flows and inflows are flipped compared to the 
travel times in a). 



Fig. 11. Velocities v ac for the average supergranule measured at the ori¬ 
gin at different solar latitudes from MDI TD (f modes and pi modes). 
Solid lines are for the average supergranule outflow region, dashed lines 
for the average supergranule inflow region. The errorbars have been 
computed from dividing the 177 datasets into eight parts. 


horizontal divergence div/* (the FWHM is about 10 Mm in both 
outflow and inflow regions). 

Differences in the vorticity strength between regions of di¬ 
vergent and convergent flows have also been reported by other 
authors who studied the statistics of vortices in solar convec¬ 
tion. Wa ng et aL| ([l995) found, on granular scales, the root mean 


square of oj z to be slightly higher in inflow regions. Potzi & 


|Brandt[ |2007| observed that vortices are strongly connected to 
sinks at mesogranular scales. Concentration of fluid vorticity 
in inflows has also been found in simulations of solar convec¬ 
tion (e.g., [Stein & Nordlund|1998j on granulation scale). How¬ 
ever, the authors did not find any preferred sign of co z . The in¬ 
creased vorticity strength in inflows might be a manifestation of 
the “bathtub effect” ( jNordiund 1985| . In that scenario, initially 
weak vorticity becomes amplified in inflows due to angular mo¬ 
mentum conservation. In the downflows that are associated with 
the horizontal inflows because of mass conservation, the vor¬ 
tex diameter is reduced since the density rapidly increases with 
depth. This further enhances the vorticity. 


rized in Fig. [12] The arrows show the horizontal velocity magni¬ 
tudes and directions at 40° latitude. The flows are dominated by 
the radial velocity component. For direct comparison, the filled 
contours give the vertical vorticity co z of the flows. In the aver¬ 
age outflow region (Fig.[T2fc), the vorticity shows a broad plateau 
region (FWHM about 13 Mm). The region of negative vorticity 
is surrounded by a ring of positive vorticity with a diameter of 
about 30 Mm. 

In contrast, the vorticity in the average inflow region 
(Fig. [12 Ip ) falls off rapidly from its narrow center (FWHM 
8 Mm). We note that the FWHM of the vorticity peak is smaller 
than for the divergence peak (about 10 Mm) but still larger 
than the FWHM of the LCT correlation measurements (roughly 
3 Mm). The peak vorticity magnitude is about twice the value of 
the outflow region (about 8 x 10 -6 s -1 anti-clockwise compared 
to 4 x 10 -6 s -1 clockwise). Like for the average outflow region, 
the central vorticity structure in the inflow region is surrounded 
by a ring of vorticity with opposite sign. The vorticity magnitude 
in the ring appears to be smaller than in Fig.[l2|i. 

These differences in the vortex structures of outflow and in¬ 
flow regions are visible at all latitudes (except at the equator, 
where we measure no net vorticity). The FWHM of the peak 
structures as well as the ratio of the peak vorticities (between 
outflow and inflow regions) are constant over the entire observed 
latitude range. Such differences do not appear in maps of the 

Article number, page 10 of 


6. Radial and tangential velocities versus radial 
distance 


Let us now look in greater detail at the isotropic part of the hor¬ 
izontal flow profile of the average supergranule. Figs. |T3|t and b 
show the azimuthal averages of v r and v t around both the aver¬ 
age supergranule outflow and inflow centers as a function of hor¬ 
izontal distance r to the outflow/inflow center at 40° latitude. In 
both cases, the magnitude of v r increases from the outflow/inflow 
center until it reaches a peak velocity (which we call v™ ax ) of 
slightly more than 300 m s -1 and -200 m s _1 , respectively at 
r = 7 Mm. The flow magnitudes then decrease and v r switches 
sign at a distance of about 14 Mm, marking the edge of the aver¬ 
age inflow/outflow region. In general, the v r curves for outflow 
and inflow regions are similar except for the difference in flow 
magnitude. 


The tangential velocity v t , on the other hand, exhibits simi¬ 
lar peak velocities for outflow and inflow regions (both « 
12 m s -1 ) but has opposite signs and reaches these peaks at 
different distances. The peak magnitude vf 1 ^ is about 26 times 
smaller than v™ ax in the outflow region and 18 times smaller in 
the inflow region. In the outflow region, the peak is located at 
r — 9 Mm, whereas it lies at r = 5 Mm around the average in¬ 
flow center. Despite the different peak locations, v t switches sign 
at a distance of about 17 Mm both around outflow and inflow 
centers. 
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Fig. 12. LCT horizontal velocity (black 
arrows) and vertical vorticity co z (filled 
contours) for the average supergranule 
at solar latitude 40°. a) Average outflow 
region, b) Average inflow region. Ar¬ 
rows for velocity magnitudes less than 
60 m s _1 are omitted. 


The different peak locations of v t possibly explain why the 
magnitude ratio of v ac between the average supergranule inflow 
and outflow region is smaller than one. The v ac measurements 
are especially sensitive to v t at r = 10 Mm (the annulus radius). 
At this distance, we have v t = 10ms _1 around outflow centers, 
but v t = -7 ms -1 around inflow centers, yielding a factor of 
-0.7. This agrees well with the ratio of the slopes for LCT v ac in 
Fig.[7]that we discussed in Sect. |4.3| 

Measuring the peak values of v r and v t at all latitudes except 
the equator leads to the following approximate relations: 


= (-0.059 ± 0.001) 


= (-0.089 ± 0.002) 


Q(/i) sin A 
Qo 

£l(A) sin A 
Qo 


for outflows, ( 6 ) 


for inflows, (7) 


where we used the differential rotation model from [Snodgrass] 
( |1984| ) to compute D(T), and Do denotes the rotation rate at 
the equator. The coefficients b out := -0.059 + 0.001 and b m := 
-0.089 ± 0.002 are remarkably constant over the whole latitude 
range from -60° to 60°, although vjf ax decreases from the equa¬ 
tor (335 m s -1 for outflows and -237 m s -1 for inflows) toward 
high latitudes (e.g., 272 m s -1 for outflows and -188 m s -1 for 
inflows at 60° north). The same trend is observed for measuring 
v™ ax west and east off disk center (e.g., 272 m s -1 for outflows 
and -187 m s -1 for inflows at 60° west), suggesting that the de¬ 
crease is a systematic center-to-limb effect. Since b out and b m 
are not affected, v™ ax is likely to suffer from the same systematic 
decrease as v™ ax . 

We can use Eqs. ^ and ([7]) to predict v / max for supergran¬ 
ules in the polar regions. Assuming b out and b m are indepen¬ 
dent of latitude even beyond A = ±60° and employing v™ ax 
from the equator where center-to-limb effects are small, the av¬ 
erage supergranule at the north pole should rotate with vf 1 ^ = 
(-20 m s -1 ) x Q(90°)/Qo for outflows and v™ x = (21 m s -1 ) x 
Q(90°)/Do for inflows. At the south pole, merely the sign of 
should change. 

Various authors have proposed models to describe vortex 
flows (e.g., |Taylor||1918| [Veronis||1959[ Simon & Weiss] 1997| ), 
introducing the (turbulent) kinematic viscosity 77 as a parame¬ 
ter that influences the tangential velocity component v t {r). In the 
Veronis model, the tangential flow is given by 

2/ 2 Q0l)sin2 

Vt 9 (A . 12 / j 2 \ ^ 

7 t 2 ( 4 + V Id 1 ) rj 


where l and d are the horizontal size and depth of the supergran¬ 
ule, respectively. This relationship is consistent with the mea¬ 
surements in Eqs. ([ 6 ]) and ([7]). Note, though, that the Veronis]( 


con¬ 
vection model does not include turbulence (beyond 77 ) or stratifi¬ 
cation. 


Taylor| ( |1918| ) presented a simple model that describes the 
decay of a narrow isolated vortex with v r = 0 due to fluid vis¬ 
cosity. In this case, the tangential velocity component is 

ar 9 

v t (r) = — exp(-r 2 /4770, (9) 

nt A 


with some constant a and the “age” of the vortex t. A least- 
squares fit of Eq. ([9]) to our measured curve v t (r) describes the 
v t profile for the average inflow surprisingly well. We can use 
4r]t ^7.8 Mm from the fit to obtain a crude estimate of the turbu¬ 
lent viscosity. By identifying the vortex age t with the supergran¬ 
ule lifetime (~1 day), we get 77 ~ (1 day/0 x 180 km 2 s -1 . This is 
similar to v alu es from the literature. For example, Duvall & Gi-| 


zon 


( 2000| and Simon & Weiss ( 1997] ) obtained 77 = 250 krrr 


using helioseismology and local correlation tracking of granules, 
respectively. The order of magnitude of our estimate for 77 also 
agrees with previous measurements of the diffusion coefficient of 
small magnetic elements ( [Jafarzadeh et al.|2014[ and references 
therein). 


7. Summary 

7.1. Validation 


We have successfully measured the horizontal divergence and 
vertical vorticity of near-surface flows in the Sun using different 
techniques (TD and LCT), as well as different instruments (HMI 
and MDI). Horizontal flow velocities from LCT compare well 
with line-of-sight Dopplergrams (correlation coefficient 0.94). 
Horizontal divergence maps from TD and LCT are in excel¬ 
lent agreement for 8 -h averaging (correlation coefficient 0.96 for 
75 < JcRq < 175). Vertical vorticity measurements from TD and 
LCT are highly correlated at large spatial scales (correlation co¬ 
efficient larger than 0.7 for kR Q < 100). 


We studied the average properties of supergranules by av¬ 
eraging over 3 000 of them in latitude strips from -60° to 60°. 
The vertical vorticity maps as measured from HMI TD and HMI 
LCT for the average supergranule agree at low and mid latitudes. 
Above ±40° latitude, however, the LCT and TD results are dif¬ 
ferent, due to geometrical center-to-limb systematic errors. Af¬ 
ter correcting for these errors using measurements at the equator 
away from the central meridian (cf. |Zhao et ah||2013| ), TD and 
LCT results agree well. For MDI, the TD maps are dominated 
by systematic errors even at low latitudes. Therefore, HMI is a 
significant improvement over MDI. 
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Fig. 13. Azimuthally averaged hori¬ 
zontal velocities around the average su¬ 
pergranule outflow and inflow centers 
at solar latitude 40°. The measurements 
were obtained from LCT of granules 
in HMI intensity images, a) Horizon¬ 
tal velocities v r and v t around the av¬ 
erage supergranule outflow center. The 
thin lines mark an estimate of the vari¬ 
ability of the data as measured from di¬ 
viding the 336 datasets into eight parts. 
For v r , the 10 cr level is shown, for v t the 
3cr level, b) As a), but around the aver¬ 
age supergranule inflow center. 


7.2. Scientific results: spatial maps of vertical vorticity 

Our findings can be summarized as follows. The root mean 
square of the vertical vorticity in a map of size -180 x 180 Mm 2 
at the equator and 8 h averaging is about 15 x 10 -6 s -1 after 
low-pass filtering (power at scales kR Q < 300). 

After averaging over several thousand supergranules, the av¬ 
erage outflow and inflow regions possess a net vertical vorticity 
(except at the equator). The latitudinal dependence of the vortic¬ 
ity magnitude is consistent with the action of the Coriolis force: 
co z (A) oc Q(T) sinT/Do- I n the northern hemisphere, horizontal 
outflows are associated with clockwise motion, whereas inflows 
are associated with anti-clockwise motion. In the southern hemi¬ 
sphere, the sense of rotation is reversed. This resembles the be¬ 
havior of high and low pressure areas in the Earth’s weather sys¬ 
tem (e.g., hurricanes). 

Vortices in the average supergranular inflow regions are 
stronger and more localized than in outflow regions. For exam¬ 
ple, at 40° latitude, the vertical vorticity is 8 x 10 -6 s -1 anti¬ 
clockwise in inflows versus 4 x 10~ 6 s -1 clockwise in outflows, 
whereas the FWHM is 8 Mm versus 13 Mm. The maximum tan¬ 
gential velocity in the average vortex is about 12 m s _1 at ±40° 
latitude, which is about 26 and 18 times smaller than the maxi¬ 
mum radial flow component for outflow and inflow regions, re¬ 
spectively. 

We have demonstrated the ability of TD and FCT to char¬ 
acterize rotating convection near the solar surface. This infor¬ 
mation can be used in the future to constrain models of turbu¬ 
lent transport mechanisms in the solar convection zone (cf., e.g., 
[Rudiger et al.|2014| ). The azimuthally averaged velocity compo¬ 
nents v r and v t for supergranular outflows and inflows at various 
latitudes are available as online data. 
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Appendix A: Ridge filters 

Prior to the travel-time measurements, the wavefield that is 
present in the Dopplergrams is filtered to select single ridges 
(the f modes or the pi modes). The goal is to capture as much 
of the ridge power as possible, even if the waves are Doppler- 
shifted due to flows. At the same time, we want to prevent power 
from neighboring ridges from leaking in and select as little back¬ 
ground power as possible. 

To construct the filter, we first measure the power spectra of 
the Dopplergrams at the equator and averaged over 60 days (59 
days) of data in the case of HMI (MDI). After further azimuthal 
averaging, we identify the frequency aj mo & where the ridge max¬ 
imum is located as a function of wavenumber k. 

The filter is constructed for each k as a plateau of width 2oj 5 
centered around the ridge maximum 6L> mo de- The lower and upper 
boundaries of the plateau we call oj^ and oj c . Next to the plateau, 
we add a transition region of width u> s i ope , which consists of a 
raised cosine function that guides the filter from one to zero, 
symmetrically around aj mo d e . The lower and upper limits of the 
filter we call oj a and cjj, respectively. 

The plateau half-width oj 6 consists of the following terms 


(o 5 (k) = 


My(k) 

2 


+ Ci) v (k) + £U CO nst> 


(A.l) 


where ujr(k) is the FWHM of the ridge (measured from the aver¬ 
age power spectra), cj v (k) = akv m2iX is the Doppler shift due to a 
hypothetical flow of magnitude v max multiplied by a scale factor 
a , and uj const is a constant term of small magnitude that broadens 
the filter predominantly at small wavenumbers. 

The width of the transition region relative to the plateau 
width is 


^slope — (A.2) 

where j is a unitless factor. 

In addition, we restrict the filter to a range of wavenumbers. 
Above and below a k interval, the filters are set to zero. The 
k limits of the interval are chosen such that the ridge power is 
roughly twice the background power. Because 6L> mo de is a func¬ 
tion of wavenumber, these limits can also be expressed as fre¬ 
quencies (O m in and 6c) max . 

Table |A.l] lists the filter parameters we chose for the f-mode 
and pi-mode ridge filters that we use throughout the paper. Note 
that we use the same filters for all latitudes and longitudes. For 
the pi modes, we also list an alternative filter that we use to dis¬ 
cuss the impact of the filter details on the travel-time measure¬ 
ments (see Appendix |C. 3 1 ). 

Appendix B: Conversion of travel times into flow 
velocities 

Point-to-point travel times r dlff (ri, r 2 ) are sensitive to flows in the 
direction of ri - r 2 . If the flow structure is known, travel times 
r dlff can be predicted with the knowledge of sensitivity kernels. 
Conversely, the velocity field can be obtained from measured 
travel times by an inversion. Such inversion processes are, how¬ 
ever, delicate to solve due to the illposedness of the problem. A 
simple way to obtain rough estimates of the flow velocity while 
avoiding inversions is the multiplication of the travel times by a 
constant conversion factor. Such a conversion factor can be cal¬ 
culated by artificially adding the signature of a uniform flow of 
known magnitude and direction to Dopplergrams. The magni¬ 
tude of the measured travel time divided by the input flow speed 


Table A.l. Parameters of the ridge filters that are used for the travel¬ 
time measurements in this paper (see text for details). The lower part of 
the table gives the filter limits at kR Q = 800. The filter limits for HMI 
and MDI are equivalent. 


Parameter 


Selected ridge 



f modes 

Pi modes 

Pi modes 



(regular) 

(alternative) 

^min/^TT 

1.75 mHz 

1.90 mHz 

1.90 mHz 

^max / 2?T 

5.00 mHz 

5.40 mHz 

5.00 mHz 

^const /27T 

0.025 mHz 

0.025 mHz 

0.030 mHz 

%nax 

500 m s -1 

500 m s _1 

500 m s -1 

a 

1.0 

1.0 

2.0 

j 

1.0 

1.0 

0.6 

oj a l2n 

2.50 mHz 

3.10 mHz 

3.03 mHz 

0Jb/2n 

2.67 mHz 

3.29 mHz 

3.20 mHz 

^mode / 2u 

2.84 mHz 

3.48 mHz 

3.48 mHz 

OJc/ln 

3.01 mHz 

3.67 mHz 

3.77 mHz 

ojd! 2n 

3.18 mHz 

3.86 mHz 

3.94 mHz 


yields the conversion factor. In the following, we describe this 
process. 

First, we create data cubes 0 v (r, t) that have Doppler-shifted 
power spectra to mimic the effect of a flow v independent of posi¬ 
tion r and time t. The data cubes are based on the noise model by 
Gizon & Birch |2004| , so signatures from flows others than v are 
not present. Following the noise model, we construct in Fourier 
space 0 v (k, oj) = y/P v (k, oj)Nq,\ (k, oj). Here k is the horizontal 
wave vector, P v is a Doppler-shifted power spectrum and, at each 
(k, co), No,i are independent complex Gaussian random variables 
with zero mean and unit variance. Employing No,i ensures that 
the values 0 v (k, oj) are uncorrelated, which means that there is no 
signal from wave scattering. We use P v (k, oj) = Po(k, oj - Soj) 
based on an average power spectrum Po that was measured from 
60 days of HMI Dopplergrams (and 59 days of MDI Doppler¬ 
grams) at the solar equator. The quantity Soj = k • v is the fre¬ 
quency shift due to a background flow v = (v*, 0) that we add. 
We construct 8 h datasets 0 v (r, t) for v x in the range between 
-1 000 and 1 000 ms -1 in steps of 100 m s -1 . For each velocity 
value, we compute 10 realizations. 

As a consistency check, we apply a second method for 
adding an artificial velocity signal to the HMI Dopplergram 
datasets. This procedure consists of tracking at an offset rate. 
The tracking parameters from |Snodgrass||l984| are modified by 
a constant corresponding to a v* velocity of -100 m s -1 and 
100 m s -1 , respectively. The tracking and mapping procedure 
is as for the regular HMI observations. We produce 112 such 
datacubes for each v* value at the solar equator. 

For both methods, the 8 h datasets are ridge-filtered like the 
normally tracked Doppler observations (f modes and pi modes). 
We measure travel times r dlff in the v direction with the pairs of 
measurement points separated by 10 Mm. This distance matches 
the separation in the r ac measurements. The reference cross¬ 
covariance C ref is taken from the regularly tracked HMI (MDI) 
observations averaged over 60 days (59 days) of data at the solar 
equator. This ensures that the artificial flow signal is captured by 
the travel-time measurements. 

The resulting r dlff values averaged over maps and datasets 
are shown for HMI in Fig. |B.1| For both f and pi modes, the 
travel times from offset tracking are systematically larger than 
for the Doppler-shifted power spectra by about 10 to 15%. In 
general, the travel-time magnitudes are larger for the f modes 
than for the pi modes for the same input velocity value. The re- 
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Fig. B.l. Point-to-point travel times 
from HMI Dopplergrams with artificial 
velocity signal. The point separation is 
10 Mm in the east-west direction, a) 
f modes, b) pi modes. The blue dots 
give the travel times from Dopplergram 
series that were constructed using the 
noise model by |Gizon & Birch| ( |2004| ). 
We applied a least-squares fit with a 
polynomial of degree three to the result¬ 
ing data (black curves). The red curves 
show the linear term of the fit. For com¬ 
parison, the black filled circles show 
travel times from HMI Dopplergrams 
that were tracked at an offset rate. 


Table B.l. Coefficients h\ of the cubic polynomial defined in Eq. ( |B.l| ) 
obtained from a least-squares fit. The coefficient hi for a distance of 
10 Mm is used to convert measured travel times into flow velocities. 
For comparison, the coefficients for other distances are also given. 


Instrument 

Distance 

Value of coefficient hi [s 2 m l ] for 



f modes 

Pi modes 

HMI 

5 Mm 

-9.61 x 10‘ 2 

-4.48 x 10" 2 


10 Mm 

-1.78 x 10- 1 

-8.98 x 10- 2 


15 Mm 

-2.20 x 10- 1 

-1.39 x 10-' 


20 Mm 

-2.15 x 10- 1 

-1.71 x 10” 1 

MDI 

5 Mm 

-7.34 x W 2 

-4.66 x 10“ 2 


10 Mm 

-1.58 x 10 -1 

-8.91 x 10“ 2 


15 Mm 

-1.97 x 10- 1 

-1.12 x 10" 1 


20 Mm 

-2.43 x 10 1 

-1.58 x 10" 1 


lation between input velocity v* and output travel time r dlff is 
linear only in a limited velocity range. Whereas this range spans 
from roughly -700 m s -1 to 700 m s -1 for the pi modes, it only 
reaches from -200 ms -1 to 200 m s -1 for the f modes. For ve¬ 
locity magnitudes larger than 700 m s _1 , the measured f-mode 
travel times even decrease. However, the supergranular motions 
that we analyze reach typical velocities of ~300 m s -1 , which is 
well below that regime. 

We applied a least-squares fit to a polynomial of degree 
three to the r dlff measurements from Doppler-shifted cubes (pink 
curve): 

rfV*) = h 3 vl + h 2 v 2 x + hv x + h 0 . (B.l) 


The linear term of the polynomial is shown for HMI as the red 
curve in Fig. B.l For the actual conversion, only the linear coef¬ 
ficient h\ is used. We obtain h\ - -0.178 s 2 m -1 for the f modes 
and h\ = -0.090 s 2 m -1 for the pi modes. For comparison, the 
coefficients h\ are listed for different distances in Table |B.l] The 
table also contains the coefficients for MDI. We convert travel 
times into velocities by multiplying the travel times by l/h \. The 
velocities obtained from converting r ac maps we call v ac . 


Appendix C: Systematic errors 

Appendix C. 1: Center-to-limb systematics 

At high latitudes, the original v ac and LCT co z maps for the av¬ 
erage supergranule show strong deviations from the azimuthally 
symmetric peak-ring structures that are visible at low latitudes. 
Considering that the magnitude of r ac and oj z is much smaller 
than the magnitude of r 01 and div/* at any latitude, it is possible 
that even a small anisotropy in the divergent flow component of 


the average supergranule is picked up by the v ac and co z measure¬ 
ments and added to the signal from the tangential flow compo¬ 
nent that we want to measure. Such anisotropies can arise from 
various origins. Among them are geometrical effects that depend 
on the distance to the disk center. 

For TD measurements, the sensitivity kernels depend on the 
distance to the limb. At 60° off disk center, r dlff sensitivity ker¬ 
nels for measurements in the direction along the limb differ 
strongly from kernels for measurements in the center-to-limb di¬ 
rection (see, e.g.,|Jackiewicz et al. 2007[ for a discussion). Ad¬ 
ditionally, there is a gradient of the root mean square travel time 
in the center-to-limb direction. 

In the case of LCT, the shrinking Sun effect causes large- 
scale gradients of the horizontal velocity (of several hundred me¬ 
ters per second) pointing towards disk center (Lisle & Toomre 
2004). This effect is presumably caused by insufficient reso¬ 
lution of the granules. Although HMI intensity and Doppler 
images have a pixel size of about 350 km at disk center, the 
point spread function has a FWHM of about twice that value. In 
Dopplergrams, the hot, bright, and broad upflows in the granule 
cores cause stronger blueshifts than the redshifts from the cooler, 
darker, and narrow downflows. Due to the insufficient resolution, 
the granules appear blue-shifted as a whole. This blue shift adds 
to the blue shift of granules that move toward the observer (i.e., 
towards disk center), giving them a stronger signal in the Dopp¬ 
lergram. Lisle & Toomre argue that LCT of Dopplergrams gives 
more weight to these granules than to those granules that move 
away from the observer. However, it is not clear what causes the 
shrinking Sun effect in LCT of intensity images. Fortunately, the 
shrinking Sun effect appears to be a predominantly large-scale 
and time-independent effect, so it can easily be removed from 
LCT velocity maps by subtracting a mean image. 

Another problem is the foreshortening. Far away from the 
disk center, the granules are not resolved as well in the center-to- 
limb direction as in the perpendicular horizontal direction. This 
introduces a dependence of the measurement sensitivity on an¬ 
gle. Indeed we measure at ±60° latitude that the radial flow com¬ 
ponent v r of the average supergranule is weaker by 15 to 20% in 
the center-to-limb direction compared to the perpendicular di¬ 
rection. This corresponds to a maximum velocity difference of 
about 50 m s -1 for outflows and 30 m s -1 for inflows. At 40° 
latitude, in contrast, this difference is less than 2% (6 ms" 1 ). 

Appendix C.2: MDI instrumental systematics 

Whereas for HMI the removal of geometrical center-to-limb ef¬ 
fects results in similar v ac peak structures in the supergranule 
outflow regions in the whole latitude range from -60° to 60°, for 
MDI the peak structures appear asymmetric and distorted even 
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Fig. C.l. Circulation velocities v ac of 
the average supergranule outflow region 
at solar latitude 40° derived from HMI 
and MDI Dopplergrams (after the cor¬ 
rection for center-to-limb systematics). 
The velocity maps have been obtained 
by applying the respective conversion 
factors from Appendix [B] to the travel 
times r ac . The limits of the colorscale 
are arbitrarily set to +15 m s _1 . 


Fig. C.2. Peak v ac values for pi modes 
using different parameter combinations 
(A, n) for the average supergranule at 
solar latitude 40°. a) In the average out¬ 
flow region, b) In the average inflow re¬ 
gion. The blue symbols give the results 
for the pi ridge filter that has been used 
throughout this paper. For the results in 
black, an alternative pi ridge filter with 
slightly different parameters was used 
(see text for details). The errorbars have 
been computed as in Fig. [ 7 ] The an¬ 
nulus radii corresponding to the vari¬ 
ous combinations (A, n) are all within 
(10.0 ± 0.5) Mm. 


after the correction. An example for f-mode TD at 40° latitude is 
shown in Fig. |C.l| Even at disk center where geometrical effects 
should not play a role, there are visible systematic features (that 
do not appear for HMI, cf. Fig. [6]). This is probably due to in¬ 
strumental effects that are specific to MDI (see, e.g., Korzennik 
|et al.|20d4| , for a discussion of instrumental errors in MDI). 


Appendix C.3: Selection of filter andr ac geometry 
parameters 

We note that the v ac velocity results for TD depend on the details 
of the ridge filter as well as the geometry parameters (A, n) of the 
r ac measurements. 

To give an idea of this, we construct an alternative pi ridge 
filter with slightly different width parameters (see Appendix [X|. 
Additionally, we select four other combinations (A ,n) of r ac 
measurements that preserve the annulus radius R , so that R is 
within (10.0 ± 0.5) Mm for all the combinations (A ,n). As for 
the standard combination (A = 10 Mm,n = 6), we use four dif¬ 
ferent angles for each additional combination. 

For all these combinations and both the standard and modi¬ 
fied pi filters, we calculated v ac for the average supergranule at 
40° latitude. The resulting peak velocities are shown in Fig. |C.2| 
for both inflow and outflow regions. Note that we did not apply 
the center-to-limb correction since it only has a weak influence 
on the peak velocity magnitude at 40° latitude. 

Evidently, the modified pi filter results in systematically 
larger v ac amplitudes. The difference with respect to the standard 
filter increases with decreasing A. For A = 10 Mm and n = 6, it is 
about 10%. This is qualitatively i n line wi th Duvall & Hanasoge 
( |2013| ). Using phase-speed filters, [Duvall & H anasoge observed 
that the strength of the travel-time signal from supergranulation 
is strongly dependent on the filter width. This shows that one 


should be careful when comparing absolute velocities from TD 
and LCT. For more reliable velocity values, an inversion of r 01 
and r ac maps would be needed. 

The comparison of different combinations (A, n) for the same 
filter shows that for n = 4, 6, and 8 the v ac amplitudes are similar, 
so selecting the combination (A = 10 Mm, n = 6), as we did for 
most of this work, appears justified. Decreasing A to about 5 Mm 
changes the peak v ac values. A possible reason is that A in this 
case becomes comparable to the wavelength of the oscillations, 
so it is harder to distinguish between flows in opposite direc¬ 
tions. For small n , on the other hand, the measurement geometry 
deviates strongly from a circular contour. This might explain the 
deviations in v ac for n = 3. 







































